home *** CD-ROM | disk | FTP | other *** search
/ MacTech 1 to 12 / MacTech-vol-1-12.toast / Source / MacTech® Magazine / Volume 06 - 1990 / 06.11 Nov 90 / NeuralNet Estimator⁄EDIT / Compute SS.c next >
Encoding:
C/C++ Source or Header  |  1989-08-03  |  6.7 KB  |  256 lines  |  [TEXT/KAHL]

  1. #include "Neural Network.h"
  2. #include <math.h>
  3.  
  4. extern FILE * Jac;
  5.  
  6. extern NeuralNet *     theNet;
  7. extern DTypeVector     yData;
  8. extern DTypeMatrix    XData;
  9. extern DTypeVector     Alpha[];
  10. extern DTypeMatrix  Jac_T;
  11. extern DTypeVector     Resid;
  12. extern DTypeVector     dSquash[];
  13. extern DTypeMatrix  Phi, T2;
  14.  
  15. /***************************************************************************/
  16. /*-------------------
  17.     Logistic squashing function and its derivative.
  18. */
  19. static        logisticSquash_dSquash(vector,dvector)
  20.     DTypeVector * vector;
  21.     DTypeVector * dvector;
  22. {
  23.     register int        i;
  24.     register DataType    temp;
  25.     register DataType     * cell;
  26.     register DataType     * dcell;
  27.     
  28.     cell = *vector->cells;
  29.     dcell = *dvector->cells;
  30.     
  31.     for(i=0; i<vector->rows; i++, cell++, dcell++)
  32.     {    temp = exp(-*cell);
  33.         *cell = 1/(1+temp);
  34.         *dcell = (*cell)*(1.0-*cell);
  35.     }
  36. }
  37.  
  38. /*-------------------
  39.     Compute output for the neural net where the input values are specified by
  40.     the Alpha[0] vector.  Also computes derivative of squash function.
  41. */
  42. static DataType        ComputeOutputJac()
  43. {
  44.     int i;
  45.     
  46.     for(i=0; i < theNet->OutLayer-1; i++)    /* skip the output layer since don't want squash */
  47.     {    Matrix_by_Vec(&theNet->W[i],&Alpha[i],&Alpha[i+1]);
  48.         logisticSquash_dSquash(&Alpha[i+1], &dSquash[i+1]);        /* logistic squash    */
  49.  
  50.     }
  51.     Matrix_by_Vec(&theNet->W[theNet->OutLayer-1], &Alpha[theNet->OutLayer-1], &Alpha[theNet->OutLayer]);
  52.         /* no squash for output layer */
  53.     return(**Alpha[theNet->OutLayer].cells);
  54. }
  55.  
  56. /*-------------------
  57.     Compute Jacobian for the SS problem where the input values are specified by
  58.     the Alpha[0] vector.  Expects activation values for each layer in 
  59.     Alpha[i], i=1,..,OutLayer
  60. */
  61. static         ComputeJacobian(jcell)
  62.     DataType    * jcell;    /* points to last element of ith column of Jacobian matrix */
  63. {
  64.     int i,j,k;
  65.     int OL;
  66.     DataType * acell;    /* pointer to elements of Alpha vector */
  67.     DataType * pcell;    /* pointer to elements of Phi matrix */
  68.     
  69.     OL = theNet->OutLayer;
  70.  
  71. /**** First, get the alpha values just prior to the output layer ****/
  72.     acell = *Alpha[OL-1].cells + Alpha[OL-1].rows-1;
  73.     for(i=0; i<Alpha[OL-1].rows; i++, acell--, jcell -= Jac_T.cols)
  74.         *jcell = *acell;    
  75.  
  76. /**** Second, calculate the Phi_OL-1 matrix and use with the second to last layer ****/
  77.     Phi.cols = theNet->W[OL-1].cols;
  78.     Matrix_by_Diag(&theNet->W[OL-1], &dSquash[OL-1], &Phi);
  79.     acell = *Alpha[OL-2].cells + Alpha[OL-2].rows-1;    /* points to the last element of the Alpha[OL-2] vector */
  80.     for(j=0; j<Alpha[OL-2].rows; j++, acell--)
  81.     {    pcell = *Phi.cells + (Phi.cols-1);    /* points to last element of Phi matrix, which is a row vector */
  82.         for(k=0; k<Phi.cols; k++, pcell--, jcell -= Jac_T.cols)
  83.             *jcell = (*acell)*(*pcell);
  84.     }
  85.  
  86. /**** Third, calculate Jacobian values for any other layers ****/    
  87.     for(i=OL-3; i>-1; i--)        /* indexing specifies the Alpha vector */
  88.     {    T2.cols = theNet->W[i+1].cols;
  89.         Matrix_by_Matrix(&Phi, &theNet->W[i+1], &T2);
  90.         Phi.cols = T2.cols;
  91.         Matrix_by_Diag(&T2, &dSquash[i+1], &Phi);
  92.         acell = *Alpha[i].cells + Alpha[i].rows-1;    /* points to the last element of the Alpha[i] vector */
  93.         for(j=0; j<Alpha[i].rows; j++, acell--)
  94.         {    pcell = *Phi.cells + (Phi.cols-1);    /* points to last element of Phi matrix, which is a row vector */
  95.             for(k=0; k<Phi.cols; k++, pcell--, jcell -= Jac_T.cols)
  96.                 *jcell = (*acell)*(*pcell);
  97.         }
  98.     }
  99. }
  100.  
  101. /*-------------------
  102.     Compute the sum of squares and Jacobian for net.
  103.     Output layer must be a singleton, so that dependent values can be
  104.     represented as a vector of scalars.  Because the SS function is only appropriate
  105.     for scalars, to use vector output variables would need to formulate a multidimensionl
  106.     sum of squares criteria.
  107. */
  108. double
  109. Compute_SS_Jac()
  110. {
  111.     int    i;
  112.     DataType * y;        /* pointer to actual dependent value */
  113.     DataType * r;        /* pointer to cells of residual vector */
  114.     DataType * jcell;    /* pointer to cells of Jacobian matrix, stored as transpose */
  115.     double SS;
  116.     DataType * handle;    /* used for pseudo vector */
  117.     
  118.     HLock(Phi.cells);
  119.     HLock(T2.cells);
  120.     HLock_Alpha_dSquash();
  121.  
  122.     Alpha[0].cells = & handle;    /* setup pseudo handle for use by Alpha[0].cells */
  123.                 /*    the idea here is to set up a pseudo vector which is a row in the input
  124.                     data matrix and then use this vector as the input vector for the 
  125.                     ComputeOutputJac function. 
  126.                 */
  127.     
  128.     y = *yData.cells;
  129.     r = *Resid.cells;
  130.     
  131.     jcell = *Jac_T.cells + (Jac_T.rows-1)*Jac_T.cols;    
  132.             /* points to first element of last row of Jac_T, that is, the last element 
  133.                 of first column of Jacobian matrix.  Since Jac_T is transpose of
  134.                 Jacobian this means it points to cell for derivative with respect
  135.                 to the last parameter for first observation.  By incrementing jcell++
  136.                 at each iteration the pointer specifies the last element in each 
  137.                 observation's column.
  138.             */
  139.     SS = 0.0;
  140.     for(i=0; i < XData.rows; i++, r++, jcell++)
  141.     {
  142.         handle = *XData.cells + i*XData.cols;
  143.         *r = y[i] - ComputeOutputJac();
  144.         SS += (double)pow( *r , 2);
  145.         ComputeJacobian(jcell);
  146.     }
  147.  
  148.     HUnlock(Phi.cells);
  149.     HUnlock(T2.cells);
  150.     HUnlock_Alpha_dSquash();
  151.     return(SS);
  152. }
  153.  
  154. /*-------------------
  155.     Logistic squashing function and its derivative.
  156. */
  157. static         logisticSquash(vector)
  158.     DTypeVector * vector;
  159. {
  160.     register int        i;
  161.     register DataType     * cell;
  162.     
  163.     cell = *vector->cells;
  164.     
  165.     for(i=0; i<vector->rows; i++, cell++)
  166.         *cell = 1/(1+exp(-*cell));
  167. }
  168.  
  169. /*-------------------
  170.     Compute output for the neural net where the input values are specified by
  171.     the Alpha[0] vector.
  172. */
  173. static DataType        ComputeOutputOnly()
  174. {
  175.     int i;
  176.     
  177.     for(i=0; i < theNet->OutLayer-1; i++)    /* skip the output layer since don't want squash */
  178.     {    Matrix_by_Vec(&theNet->W[i],&Alpha[i],&Alpha[i+1]);
  179.         logisticSquash(&Alpha[i+1]);        /* logistic squash    */
  180.     }
  181.     Matrix_by_Vec(&theNet->W[theNet->OutLayer-1], &Alpha[theNet->OutLayer-1], &Alpha[theNet->OutLayer]);
  182.         /* no squash for output layer */
  183.     return(**Alpha[theNet->OutLayer].cells);
  184. }
  185.  
  186. /*-------------------
  187.     Compute the sum of squares
  188. */
  189. double
  190. Compute_SS()
  191. {
  192.     int    i;
  193.     DataType * y;        /* pointer to actual dependent value */
  194.     double SS;
  195.     DataType * handle;    /* used for pseudo vector */
  196.     double r = 0.0;
  197.  
  198.     HLock_Alpha();
  199.     
  200.     Alpha[0].cells = & handle;
  201.     y = *yData.cells;
  202.     SS = 0.0;
  203.     for(i=0; i < XData.rows; i++)
  204.     {
  205.         handle = *XData.cells + i*XData.cols;
  206.         r = y[i] - ComputeOutputOnly();
  207.         SS += (double)pow( r , 2);
  208.     }
  209.  
  210.     HUnlock_Alpha();
  211.     
  212.     return(SS);
  213. }
  214.  
  215. /*-------------------
  216. */
  217. static         HLock_Alpha_dSquash()
  218. {
  219.     int i;
  220.     
  221.     for(i=0; i<theNet->OutLayer; i++)
  222.     {    HLock(Alpha[i].cells);
  223.         HLock(dSquash[i].cells);
  224.     }
  225. }
  226.  
  227. static         HUnlock_Alpha_dSquash()
  228. {
  229.     int i;
  230.     
  231.     for(i=0; i<theNet->OutLayer; i++)
  232.     {    HUnlock(Alpha[i].cells);
  233.         HUnlock(dSquash[i].cells);
  234.     }
  235. }
  236.  
  237. static         HLock_Alpha()
  238. {
  239.     int i;
  240.     
  241.     for(i=0; i<theNet->OutLayer; i++)
  242.     {    HLock(Alpha[i].cells);
  243.     }
  244. }
  245.  
  246. static         HUnlock_Alpha()
  247. {
  248.     int i;
  249.     
  250.     for(i=0; i<theNet->OutLayer; i++)
  251.     {    HUnlock(Alpha[i].cells);
  252.     }
  253. }
  254.  
  255.  
  256.